El archivo manchasolares.txt, contiene el número observado de manchas solares en el Sol en cada mes desde enero de 1749. El archivo contiene dos columnas de números, la primera es el mes y la segunda el número de manchas solares.
(a) Escribe un programa que lea los datos y haga una gráfica de las manchas solares en función del tiempo.
from pylab import*
import matplotlib.pyplot as plt
import numpy as np
#Me debo acostumbrar a esto, porque si luego hay cosas que se llaman iguales, se confunde.
#Cuando estemos en duda preguntamos a Google
file = open('manchasolares.txt', 'r')
#Lo que hacemos es abrir el archivo y ponerlo en read mode.
#Lit es una analogia a la vida real, tomas un libro, lo abres, lo lees y despues tu decides que hacer con el
#escribilre (modificarlo, varias cosas), y despues de abrir el libro tienes que cerrarlo.
line = file.readline()
# Aquí la instruccion se la asignamos a la variable line, es decir, va a leer las lineas del archivo
#Nosotros sabemos que las lineas son del tipo # #, donde gato son cualquier número, y son separados por el
#caracter \t (un doble espacio quiero creer yo), eso es lo que me sugirio el copilot (no es como que yo supiera que era ese)
#*cries in regex*
Meses = []
Manchas = []
#Lo tipico, debo crear unas listas vacias, para que pueda distinguir el tipo de datos que voy guardando
for line in file:
items = line.split('\t')
Meses.append(float(items[0]))
Manchas.append(float(items[1]))
file.close()
figure(figsize=(10, 10), dpi=800)
#No me acordaba de esto, pero es para hacer la figura más grande
plt.plot(Meses, Manchas)
plt.grid(True)
plt.xlabel('Meses')
plt.ylabel('Manchas')
plt.title('Manchas solares en función del tiempo')
plt.show()
(b) Modifica tu programa para mostrar solo los primeros 1000 datos (experimentales) en la gráfica.
file = open('manchasolares.txt', 'r')
line = file.readline()
Meses = []
Manchas = []
for line in file:
items = line.split('\t')
Meses.append(float(items[0]))
Manchas.append(float(items[1]))
file.close()
figure(figsize=(10, 10), dpi=800)
plt.plot(Meses[:1000], Manchas[:1000],color = 'red') #lit, solo fue ponerle tope a que datos deberia considerar
plt.grid(True)
plt.xlabel('Meses')
plt.ylabel('Manchas')
plt.title('Manchas solares en función del tiempo (1000 datos)')
plt.show()
c) Modifica nuevamente tu programa para calcular y graficar la media (promedio) móvil de los datos, definida por:
$$Y_k = {1\over2r} \sum_{m=-r}^r y_{k+m}$$donde $r = 5$ (en este caso) y $y_k$ son los números de manchas solares.
El programa debe graficar tanto los datos originales como la media móvil en el mismo gráfico, sólo sobre los primeros 1000 datos.
file = open('manchasolares.txt', 'r')
line = file.readline()
Meses = []
Manchas = []
for line in file:
items = line.split('\t')
Meses.append(float(items[0]))
Manchas.append(float(items[1]))
file.close()
#Okay, mantenemos esta parte intacta, lo unico que hace es darnos los datos.
#Definimos nuevas cosas
r = 5
MediaMovil = []
for k in range(0,1000):
#1000 datos asi es
Suma = 0
for m in range(-r,r+1):
if k+m >= 0 and k+m <= 1000:
#ocupamos una condicion, pues al solo trabajar con 1000 datos,
#no tenemos acceso a el dato -5,-4-3 y asi, o el 1001, 1002, 1003,...
Suma = Suma + Manchas[k+m]
MediaMovil.append(Suma/(2*r)) #podriamos poner solo 1/10, pero asi es mas general
#al terminar el for para m, pues tenemos nuestro Y_k
#Ploteamos asi es
figure(figsize=(10, 10), dpi=800)
plt.plot(Meses[:1000], Manchas[:1000],'green',label = 'Manchas solares')
plt.plot(Meses[:1000], MediaMovil[:1000],color = 'orange',label = 'Media movil')
plt.grid(True)
plt.xlabel('Meses')
plt.ylabel('Manchas/Media Movil')
plt.title('Manchas solares en función del tiempo (1000 datos)')
plt.legend(loc="center right")
plt.show()
Aunque la función plot está diseñada principalmente para hacer gráficos $xy$ estándar, se puede adaptar para otros tipos de gráficas también.
a) Grafica la llamada curva deltoide, definida paramétricamente por las ecuaciones:
$$x = 2 cos θ + cos 2θ, \; \; \; \; \; \; y = 2 sin θ − sin 2θ$$donde $0\le\theta<2\pi$. Toma un conjunto de valores para $\theta$ entre $0$ y $2\pi$ y calcula $x(\theta)$ e $y(\theta)$ usando las ecuaciones anteriores, para posteriormente graficar $y$ como función de $x$.
#The easy and lazy way
theta = linspace(0,2*np.pi,1000)
x = 2*np.cos(theta)+np.cos(2*theta)
y = 2*np.sin(theta)-np.sin(2*theta)
figure(figsize=(10, 10), dpi=800)
plt.plot(x,y)
plt.grid(True)
plt.title(r'x($\theta$) vs y($\theta$)')
plt.xlabel(r'x($\theta$)')
plt.ylabel(r'y($\theta$)')
show()
b) Usando este mismo enfoque, se puede hacer una gráfica polar $r=f(\theta)$ para alguna función $f$ calculando $r$ para un rango de valores de $\theta$ y luego convirtiendo $r$ y $\theta$ a coordenadas cartesianas usando las ecuaciones estándar:
$$x = r\cos\theta, \; \; \; \; y = r\sin\theta $$Utiliza este método para trazar la espiral Galileana: $r = \theta^2$ para $0\le\theta\le10\pi$.
#Bueno, lo que vamos a recibir va a ser una función en terminos de theta
#entonces, lo único que vamos a usar sera una r que en si tendriamos que substituir en
# x = r*cos(theta) y y = r*sin(theta)
def OsoPolar(r,start,finish,puntos,Titulo):
theta = linspace(start,finish,puntos)
x = r*cos(theta)
y = r*sin(theta)
plot(x,y)
grid(True)
title(Titulo)
xlabel(r'x($\theta$)')
ylabel(r'y($\theta$)')
show()
r_b = theta**2
OsoPolar(r_b,0,10*np.pi,1000,'Espiral Galileana')
c) Con el mismo método, haz una gráfica polar de la "función de Fey'':
$$r = e^{\cos\theta} - 2 \cos 4\theta + \sin^5 \frac{\theta}{12}$$en el rango $0\le\theta\le24\pi$.
r_c = np.exp(np.cos(theta))-2*np.cos(4*theta)+(np.sin(theta/12))**5
OsoPolar(r_c,0,24*np.pi,1000,'Función de Fey')
En un especial halloween de los Simpsons, Homero viaja al pasado a la era jurásica y por un pequeño descuido cambia el curso de toda la historia en el futuro. Lo anterior es una referencia al cuento A Sound of Thunder de 1952 del escritor de ciencia ficción Ray Bradbury en donde lo que se destruye en el pasado es una mariposa.
Esta idea fue retomada y popularizada por el físico Edward Lorenz en 1972 cuando dio una conferencia en la Asociación Estadounidense para el Avance de la Ciencia (American Association for the Advancement of Science) titulada "¿El aleteo de una mariposa en Brasil puede provocar un tornado en Texas?". Dando paso al concepto del efecto mariposa, del cual probablemente hayas oído hablar y que es el ejemplo clásico de caos determinista en sistemas climáticos.
El caos determinista también aparece en muchos sistemas físicos más complejos, incluyendo especialmente la dinámica de fluidos. Debido a su naturaleza aparentemente aleatoria, el comportamiento de los sistemas caóticos es difícil de predecir y se ve fuertemente afectado por pequeñas perturbaciones en las condiciones iniciales. Uno de los ejemplos más famosos del fenómeno del caos determinista es sin duda el mapeo logístico, que es un sistema matemático muy simple, definido por la ecuación:
$$x_{n+1} = rx_n(1-x_n)$$$\;$
Para un valor dado de la constante $r$, se toma un valor de $x_n$ (digamos $\frac{1}{2}$) y se introduce en el lado derecho de esta ecuación y regresa un valor de $x_{n+1}$. Luego se toma ese valor y se vuelve a introducir en el lado derecho, lo que da otro valor, y así sucesivamente. Esto es un mapeo iterativo. Se continua haciendo la misma operación una y otra vez sobre su valor de $x_n$, y entonces sucede una de las tres siguiente situaciones:
Responde a las siguientes preguntas:
a) Apoyate en el programa que vimos en clase y escribe un programa que muestre el comportamiento del mapeo logı́stico mediante una gráfica.
def Logistico(r,x):
return r*x*(1-x)
def Mapeo_Logistico(r, x0, n):
dominio = []
Logis = [x0]
for i in range(n):
x_n = r * Logis[i] * (1 - Logis[i]) #nuestro x[n+1]
Logis.append(x_n)
dominio.append(i)
return Logis, dominio
def Mapeo_Logistico_Ploteada(r, x0, n):
dominio = []
Logis = [x0]
Logis[0] = x0
for i in range(n):
x_n = r * Logis[i] * (1 - Logis[i]) #nuestro x[n+1]
Logis.append(x_n)
dominio.append(i)
dominio.append(n)
plt.plot(dominio, Logis, '.')
plt.xlabel('n')
plt.ylabel('x[n]')
plt.title(f'Mapeo Logistico con r={r} y x0={x0}')
plt.show()
Mapeo_Logistico_Ploteada(3.59,0.2,1000)
Mapeo_Logistico_Ploteada(1,0.2,1000)
b) De acuerdo a tu gráfica, ¿a qué valor de $r$ el sistema pasa de un comportamiento ordenado (puntos fijos u orbitas) a un comportamiento caótico? A este punto a veces se le llama "el borde del caos".
Hint: Esto es lo que debes hacer para hacer tu programa:
Para un valor dado de $r$, comienza con $x_0=\frac{1}{2}$ e itera la ecuación del mapeo logístico mil veces. Eso le dará la oportunidad de establecerse en un punto fijo o en una orbita de algún periodo. Luego ejecuta otras mil iteraciones y grafica los puntos $(r,x_{\infty})$ en una gráfica donde el eje horizontal es $r$ y el eje vertical es $x_{\infty}$. Puedes usar la función plot con las opciones "ko" o "k." para dibujar una gráfica con puntos, uno para cada punto, o puedes usar la función scatter para dibujar un diagrama de dispersión (que siempre usa puntos). Repite todo el cálculo para valores de $r$ de 1 a 4 en pasos de 0.01, graficando los puntos para todos los valores de $r$ en la misma figura. Tu programa debería generar la distintiva gráfica que parece un árbol inclinado hacia un lado. Esta famosa imagen se llama Gráfica de Feigenbaum[1], en honor a su descubridor Mitchell Feigenbaum.
def ML2(r, x0, n):
Logis = [x0]
for i in range(n):
x_n = r * Logis[i] * (1 - Logis[i])
Logis.append(x_n)
return Logis
r_bolsadevalores = np.linspace(1, 4, 1000) #de 0.001 en 0.001
x0 = 0.5
erre = []
x_infinito = []
for r in r_bolsadevalores:
equis = ML2(r, x0, 1000)
erre.append([r] * 1000)
x_infinito.append(equis[-1000:])
erre = np.array(erre).flatten() #nos ahorramos un for split append xD
x_infinito = np.array(x_infinito).flatten()
plt.figure(figsize=(10, 6))
plt.plot(erre, x_infinito, '.k', alpha=0.25, markersize=0.5) #._. ko
plt.xlabel('Valor de r')
plt.ylabel('Valor de x∞')
plt.title('Árbol de Feigenbaum')
plt.xlim(1, 4)
plt.ylim(0, 1)
plt.show()
c) Opcional (para 1.5 puntos extra): Hay otra forma para calcular el diagrama de Feigenbaum, que puede ser más claro y rápido, dado que hace uso de la capacidad de Python para realizar aritmética con arreglos completos.
Crea un arreglo r que contenga cada valor distinto de $r$, [1.0, 1.01, 1.02, ... ]. Después, crea otro arreglo x del mismo tamaño para guardar los valores correspondientes de $x$, establecidos inicialmente en $0.5$; finalmente realiza una iteración del mapeo logístico para todos los valores de $r$ a la vez, con una sola instrucción de la forma x = r*x*(1-x) y comparala con tu programa anterior[2].
#._.
1. ^ A veces también se suele llamar gráfica de higuera, debido al juego de palabras que surge del hecho de que parece un árbol y Feigenbaum significa "higuera" en alemán.
2. ^ Debido a la velocidad con la que Python puede realizar cálculos con arreglos, este método debería ser significativamente más rápido que el método más básico anterior.
El conjunto de Mandelbrot, llamado así por su descubridor, el matemático francés Benoit Mandelbrot, es un fractal; un objeto matemático infinitamente ramificado que contiene estructura dentro de estructura dentro de estructura, tan profundamente como queramos mirar. La definición del conjunto de Mandelbrot en términos de números complejos es la siguiente.
Consideremos la ecuación: $$z_{n+1} = z_n^2 + c, \; \; \; \; \; \; \; \; \; z_n, c \in \mathbb{C}$$ $\;$
donde $z_n$ es un número complejo y $c$ es una constante compleja. De manera, muy similar al mapeo logístico, la definición del conjunto de Mandelbrot implica la iteración repetida de esta ecuación para cualquier valor dado de $c$. La ecuación convierte un número de entrada $z_n$ en un número de salida $z_{n+1}$; de tal manera que se toma un valor inicial de $z_0$ y se introduce en la ecuación para obtener un nuevo valor $z_1$, luego tomamos ese valor y lo volvemos a introducir para obtener $z_2$ y así sucesivamente.
Así, el conjunto de Mandelbrot es el conjunto de puntos del plano complejo que satisface la siguiente definición:
_Para un valor dado de $c$ y la condición inicial $z_0=0$; si al iterar repetidamente la ecuación, la magnitud del valor resultante es mayor a dos (i.e. $z_{\infty} > 2$), entonces el punto del plano complejo para ese valor $c$ no está en el conjunto de Mandelbrot; de lo contrario, si está en el conjunto._
Para utilizar esta definición, en principio, habría que iterar infinitas veces para demostrar que un punto está en el conjunto de Mandelbrot, ya que un punto está en el conjunto sólo si la iteración nunca pasa de $z_n=2$. Sin embargo, en la práctica, simplemente se realiza una gran cantidad de iteraciones (digamos 100) y si $z_n$ no ha excedido 2 en ese momento, entonces lo consideramos suficientemente bueno.
a) Escribe un programa para crear una imagen del conjunto de Mandelbrot.
Realiza la iteración para todos los valores de $c=x+iy$ en una cuadrícula de $N \times N$ que abarque la región donde $-2\le x\le 2$ y $-2\le y\le 2$. Haz una gráfica de densidad (density plot) en el que los puntos de la cuadrícula dentro del conjunto de Mandelbrot estén coloreados en negro y los de afuera estén coloreados en blanco.
#in progress...
Sugerencia: Probablemente te resulte útil comenzar con una cuadrícula muy simple, es decir, con un valor pequeño de $N$ (quizás $N=100$) para que tu programa se ejecute rápidamente mientras lo pruebas. Una vez que te asegures de que funciona correctamente, aumenta el valor de $N$ para producir una imagen final de alta calidad de la forma del conjunto.
Si te aburrió lo anterior (o te resultó demasiado fácil), puedes programar otra variante del mismo ejercicio que puede producir imágenes sorprendentes. En lugar de colorear los puntos solo en blanco o negro, colorea los puntos de acuerdo con el número de iteraciones de la ecuación antes de que $z_n$ sea mayor que 2 (o bien el número máximo de iteraciones si es que $z_n$ nunca llega a ser mayor que 2). Si usas alguno de los esquemas más coloridos que Python proporciona para las gráficas de densidad, como "hot" o "jet", puedes crear algunas imágenes muy espectaculares. Otra variante interesante es colorear según el logaritmo del número de iteraciones, lo que ayuda a revelar parte de la estructura más fina fuera del conjunto.
#in progress...
Durante el curso, hemos escrito un par de programas para calcular e imprimir el factorial de un número ingresado por el usuario, usando valores enteros y no de punto flotante. Usa alguno de esos programas para calcular el factorial de 200 y modificalo para usar variables de punto flotante y calcula nuevamente el factorial de 200.
Qué es lo que encuentras?, explica.
def factorial_flotante(n):
if n == 0:
return 1.0
else:
return float(n) * factorial_flotante(float(n - 1))
def factorial(n):
if n == 0:
return 1
else:
return n * factorial(n - 1)
print(factorial(200))
print(factorial_flotante(200))
788657867364790503552363213932185062295135977687173263294742533244359449963403342920304284011984623904177212138919638830257642790242637105061926624952829931113462857270763317237396988943922445621451664240254033291864131227428294853277524242407573903240321257405579568660226031904170324062351700858796178922222789623703897374720000000000000000000000000000000000000000000000000 inf
Brutal, pues como los float tienen más cifras, se descontrola el asunto y en realidad, no es porque tengan más cifras sino que tienen un límite estos amigazos no como los enteros. Overflow, pero porque ya no dan para más, tienen un tope y python nos lo pone como infinito.
Considera una ecuación cuadrática $ax^2+bx+c=0$ que tiene soluciones reales.
a) Escribe un programa que tome como entrada tres números, $a$, $b$ y $c$, e imprima las dos soluciones de la ecuación cuadrática $ax^2+bx+c=0$; usando la fórmula estándar:
$$x_{1,2} = {-b \pm \sqrt{b^2-4ac}\over2a}$$$\;$
Utiliza tu programa para calcular las soluciones de $0.001 x^2 + 1000 x + 0.001 = 0$.
def cuadratica(a,b,c):
x1 = (-b + np.sqrt(b**2 - 4*a*c))/(2*a)
x2 = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)
return x1,x2
cuadratica(0.001,1000,0.001)
(-9.999894245993346e-07, -999999.999999)
b) Existe otra forma de escribir las soluciones de una ecuación cuadrática. Demuestra que multiplicando la parte superior e inferior de la solución anterior por $-b\pm \sqrt{b^2-4ac}$, las soluciones también se pueden escribir como:
$$x_{1,2} = {2c\over-b\pm \sqrt{b^2-4ac}}$$$\;$
Agrega unas lineas a tu programa para imprimir también estos valores, además de los anteriores y entonces usa nuevamente el programa para resolver $\; 0.001 x^2 + 1000 x + 0.001= 0$. ¿Qué es lo que ves? ¿Cómo lo explicas?.
def cuadratica2(a,b,c):
x11 = (-b + np.sqrt(b**2 - 4*a*c))/(2*a)
x21 = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)
x12 = (2*c)/(-b + np.sqrt(b**2 - 4*a*c))
x22 = (2*c)/(-b - np.sqrt(b**2 - 4*a*c))
return x11,x21,x12,x22
cuadratica2(0.001,1000,0.001)
(-9.999894245993346e-07, -999999.999999, -1000010.5755125057, -1.000000000001e-06)
c) Usando lo que has aprendido, escribe un nuevo programa que calcule con precisión ambas raíces de la ecuación cuadrática en todos los casos.
def Sol_Cuadraticas(a, b, c):
discriminante = b**2 - 4*a*c
if discriminante > 0:
x1 = (-b + np.sqrt(discriminante)) / (2*a)
x2 = (-b - np.sqrt(discriminante)) / (2*a)
return x1, x2
elif discriminante == 0:
x = -b / (2*a)
return x, x
else:
Re = -b / (2*a)
Im = np.sqrt(abs(discriminante)) / (2*a)
x1 = complex(Re, Im) #no conocia esta de complex, se ve bacano
x2 = complex(Re, -Im)
return x1, x2
#Lit solo era separar por lo que puede ser el discrimante y a partir de eso, hacer las condiciones.
Sol_Cuadraticas(9,-900,99991)
((50+92.79068439833338j), (50-92.79068439833338j))
Este es un buen ejemplo de cómo las computadoras no siempre funcionan como se espera. Si simplemente se aplica la fórmula estándar para la ecuación cuadrática, la computadora a veces obtendrá una respuesta incorrecta. En la práctica, el método que has desarrollado aquí es la forma correcta de resolver una ecuación cuadrática en una computadora, aunque es más complicado que la fórmula estándar. Si estuvieras escribiendo un programa que implicara resolver muchas ecuaciones cuadráticas, este método pod ría ser un buen candidato para convertirlo en una función: podrías poner los detalles del método dentro de dicha función para ahorrarte la molestia de seguirlo paso a paso, cada vez que tengas una nueva ecuación que resolver.
Supongamos que tenemos una función $f(x)$ y queremos calcular su derivada en un punto $x_0$. Esto lo podemos hacer con lápiz y papel si conocemos la forma matemática de la función, o podemos hacerlo en la computadora haciendo uso de la definición de la derivada:
$${df \over dx} = \lim_{\delta x \to 0} {f(x_0+\delta)-f(x_0)\over \delta x}$$$\;$
En la computadora no podemos tomar el límite cuando $\delta x$ llega a cero, pero podemos obtener una aproximación razonable simplemente haciendo que $\delta x $ sea pequeño.
a) Escribe un programa que defina una función f(x) y que devuelva el valor $x(x-1)$, luego calcula la derivada de la función en el punto $x_0=1$ usando la fórmula anterior con $\delta x =10^{-2}$. Calcula analíticamente el valor real de la misma derivada y compáralo con la respuesta que da tu programa. Ambos no estarán perfectamente de acuerdo.
¿Por qué no?
def f(x):
return x*(x-1)
def derivada(f,x0,delta):
return (f(x0+delta)-f(x0))/delta
derivada(f,1,10**-2)
1.010000000000001
Pues esta cosa tiene limitantes, y usa flotantes con fin, no como nosotros que usamos enteros (pero en retrospectiva, si nosotros usarmos flotantes, el límite es la imaginación)
b) Repite el cálculo para $\delta x =10^{-4}$, $10^{-6}$, $10^{-8}$, $10^{-10}$, $10^{-12}$ y $10^{-14}$. Deberías ver que la precisión del cálculo inicialmente mejora a medida que $\delta x $ se hace más pequeña, pero luego vuelve a empeorar.
¿Por qué pasa esto?
Bueno, pues como ya habiamos dicho en la parte del facotrial que usa flotantes, son números finitos no exactos,si continuamos usando números más pequeños llega un punto en el que dejan de servir, porque la computadora no los puede hacer infinitos, tiene memoria limitada.
for i in range(4,14+2,2):
print(derivada(f,1,10**-i), i)
1.0000999999998899 4 1.0000009999177333 6 1.0000000039225287 8 1.000000082840371 10 1.0000889005833413 12 0.9992007221626509 14
Estudiaremos las derivadas numéricas con más detalle más adelante en el curso, donde examinaremos¸ técnicas para abordar estos problemas y maximizar la precisión de nuestros cálculos.